Using Basemap with Diana

While Diana has map data that is reasonable for interactive use, higher resolution data is provided by Basemap. Ideally, there would be a way to combine the map data from Basemap with the plots from Diana. This document shows one way to do this for raster images.

Creating a Basemap Plot

To begin with, we import matplotlib and the Basemap class from the basemap module. We also create a function that lets us save plots as images and embed them in an ipython notebook.



In [1]:

    
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap

from IPython.core.display import Image
def embed_fig(fig, **kwargs):
    fig.savefig("/tmp/temp.png", **kwargs)
    return Image(data=open("/tmp/temp.png").read(), format="png", embed=True)

We create a Basemap object using a number of familiar parameters to describe a projection, since this class is based on the PROJ.4 library. We use the class to obtain coordinates for the extent of the area to be plotted.



In [2]:

    
bm = Basemap(projection="ortho", lat_0 = 60, lon_0 = 15, resolution = "i")
bmx1, bmy1 = bm(-20, 35)
bmx2, bmy2 = bm(20, 65)

We set the extent to be plotted by assigning these values to attributes in the Basemap object.



In [3]:

    
bm.llcrnrx = bmx1
bm.urcrnrx = bmx2
bm.llcrnry = bmy1
bm.urcrnry = bmy2

Finally, we create a plot and draw coastlines, countries and other navigational aids before creating and embedding an image.



In [4]:

    
plot = bm.drawcoastlines()
bm.drawcountries()
bm.drawparallels(range(-90,90,10))
bm.drawmeridians(range(-180,180,10))
embed_fig(plot.figure, bbox_inches="tight", pad_inches=0)

Since we want to overlay data onto this image, we need to ensure that we can use the same projection using Diana to create an image of the same size.

Creating a Diana Plot

We import the modules needed to access Diana's functionality, plus the pyproj module that we use to describe projections, and also the QImage and QPainter classes from PyQt4.QtGui that we will use to perform some image composition.



In [5]:

    
from metno import bdiana, metlibs, plotcommands
from metno.ipython_extensions import embed
import pyproj
from PyQt4.QtGui import QImage, QPainter

As usual, we create a BDiana object and call its setup method. Omitting the file name causes a default setup file to be used.



In [6]:

    
b = bdiana.BDiana()
b.setup()

The Basemap image was saved to a file, so we load it into a QImage object and query its size. We will need to create an image of this size using Diana.



In [7]:

    
bm_image = QImage("/tmp/temp.png")
pw, ph = bm_image.width(), bm_image.height()

Using the projection from the Basemap plot, we create a projection for use with Diana and use it to determine the rectangle that describes the plot area.



In [8]:

    
p = pyproj.Proj(bm.proj4string)
x1, y1 = p(-20, 35)
x2, y2 = p(20, 65)

We create an Area object using the projection and rectangle, and a Map which we simply fill with a transparent background (the 0 value indicates an alpha component of zero).



In [9]:

    
a = plotcommands.Area(proj4string = p.srs, rectangle = [x1, x2, y1, y2])
m = plotcommands.Map(map="Gshhs-Auto", backcolour="white:0", land="off",
                     contour="off", lat="off", lon="off")
#m.setOption("cont.colour", "red")
#m.setOption("cont.linetype", "dot")

We create a Field plot command to describe the field that we wish to show, setting its alpha option to 32 (from a range of 0 to 255) to indicate that it should be partially transparent.



In [10]:

    
f = plotcommands.Field(model="HIRLAM.12KM.00", plot="air_temperature_2m",
                       plottype="contour", palettecolours="standard",
                       linetype="solid", linewidth=1, alpha=32, repeat=1)
f.setOption("line.interval", 2)
f.setOption("colour", "black")
times = b.getFieldTimes(f.options["model"], f.options["plot"])

Having obtained a list of available times for the field, we set the plot commands, set the plot time, and create an image of the same size as the Basemap image.



In [11]:

    
b.setPlotCommands([a.text(), f.text(), m.text()])
b.setPlotTime(times[-1])
diana_image = b.plotImage(pw, ph)
embed(diana_image)

Combining the Images

To combine the images, we make a copy of the QImage object that we loaded the Basemap image into and use a QPainter object to paint the Diana image over it.



In [12]:

    
combined_image = bm_image.copy()

painter = QPainter()
painter.begin(combined_image)
painter.drawImage(0, 0, diana_image)
painter.end()

Since the field plot on the Diana image is partially transparent and the background is completely transparent, the result is an image that combines the field data from the Diana image with the features from the Basemap image.



In [13]:

    
embed(combined_image)